0.5.1 SeuratObject

Single-cell data can be provided in a number of formats, one of which is the sparse matrix format. These can be read into R using the Seurat function mentioned above and supplying the path to the directory containing the files:

setty_rep1_matrix <- Read10X("setty_rep1/")

## Warning in readLines(con = barcode.loc): incomplete final line found on
## 'setty_rep1//barcodes.tsv.gz'

The expression matrix is loaded as a dgCMatrix object, which is a compressed matrix format that saves space by only reporting non-zero elements:

class(setty_rep1_matrix)

## [1] "dgCMatrix"
## attr(,"package")
## [1] "Matrix"

head(setty_rep1_matrix[c(1:10),c(1:10)])

## 6 x 10 sparse Matrix of class "dgCMatrix"

##   [[ suppressing 10 column names 'AAACCTGAGAAACGCC', 'AAACCTGAGCAATCTC', 'AAACCTGCAAGCGAGT' ... ]]

##                                
## DDX11L2     . . . . . . . . . .
## DDX11L1     . . . . . . . . . .
## WASH7P      . . . . . . . . . .
## MIR6859-1   . . . . . . . . . .
## MIR1302-2HG . . . . . . . . . .
## MIR1302-2   . . . . . . . . . .

Note that the matrix is still formatted as usual: with columns being cell IDs (typically encoded as cell barcode sequences) and rows representing genes:

head(colnames(setty_rep1_matrix))

## [1] "AAACCTGAGAAACGCC" "AAACCTGAGCAATCTC" "AAACCTGCAAGCGAGT" "AAACCTGCACTTAAGC"
## [5] "AAACCTGGTAAGTAGT" "AAACCTGGTACCCAAT"

head(rownames(setty_rep1_matrix))

## [1] "DDX11L2"     "DDX11L1"     "WASH7P"      "MIR6859-1"   "MIR1302-2HG"
## [6] "MIR1302-2"

We can use this matrix to create a SeuratObject to store expression as well as other metrics and metadata. We will use the CreateSeuratObject() function, setting min.cells = 1 to only include genes that are expressed in at least one cell. You may read more about this function accessing the help page (?CreateSeuratObject).

setty_seurat <- CreateSeuratObject(setty_rep1_matrix,
                                   project = "healthyBM",
                                   min.cells = 1)

## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')

setty_seurat

## An object of class Seurat 
## 35102 features across 6381 samples within 1 assay 
## Active assay: RNA (35102 features, 0 variable features)

Observe that the dataset contains 39306 genes (i.e. features) and 6998 cells (i.e. samples). The SeuratObject only contains one assay, named “RNA”. If more modalities or other data types were supplied, they would be stored as additional assays in the same SeuratObject.

The SeuratObject includes multiple slots where different types of information can be stored. The @ operator opens up the different slots in the object, including the assays slot, where the “RNA” assay is contained. Inside it, there other slots featuring assay-specific information. Here’s how the hierarchy works:

setty_seurat@assays

## $RNA
## Assay data with 35102 features for 6381 cells
## First 10 features:
##  ENSG00000238009, CICP27, ENSG00000268903, ENSG00000241860,
## ENSG00000228463, ENSG00000290385, ENSG00000230021, ENSG00000235146,
## MTND1P23, MTND2P28

setty_seurat@assays$RNA

## Assay data with 35102 features for 6381 cells
## First 10 features:
##  ENSG00000238009, CICP27, ENSG00000268903, ENSG00000241860,
## ENSG00000228463, ENSG00000290385, ENSG00000230021, ENSG00000235146,
## MTND1P23, MTND2P28

head(setty_seurat@assays$RNA@counts[c(1:10),c(1:10)])

## 6 x 10 sparse Matrix of class "dgCMatrix"

##   [[ suppressing 10 column names 'AAACCTGAGAAACGCC', 'AAACCTGAGCAATCTC', 'AAACCTGCAAGCGAGT' ... ]]

##                                    
## ENSG00000238009 . . . . . . . . . .
## CICP27          . . . . . . . . . .
## ENSG00000268903 . . . . . . . . . .
## ENSG00000241860 . . . . . . . . . .
## ENSG00000228463 . . . . . . . . . .
## ENSG00000290385 . . . . . . . . . .

The object also contains cell-level information, known as metadata, which can be accessed as follows:

head(setty_seurat[[]])

##                  orig.ident nCount_RNA nFeature_RNA
## AAACCTGAGAAACGCC  healthyBM       6164         2808
## AAACCTGAGCAATCTC  healthyBM       4550         2258
## AAACCTGCAAGCGAGT  healthyBM       8697         3406
## AAACCTGCACTTAAGC  healthyBM       4325         1946
## AAACCTGGTAAGTAGT  healthyBM       2144         1370
## AAACCTGGTACCCAAT  healthyBM      10389         3940

These statistics are calculated by default upon SeuratObject creation, but more columns will be added to this table as we perform additional analyses.

0.5.2 SingleCellExperiment

The loaded dgCMatrix object can be made into a SingleCellExperiment object using the SingleCellExperiment() function. In this case, we will supply a list object in which every element will be one assay. This allows us to set assay names:

setty_sce <- SingleCellExperiment(list(RNA=setty_rep1_matrix))

class(setty_sce)

## [1] "SingleCellExperiment"
## attr(,"package")
## [1] "SingleCellExperiment"

setty_sce

## class: SingleCellExperiment 
## dim: 62754 6381 
## metadata(0):
## assays(1): RNA
## rownames(62754): DDX11L2 DDX11L1 ... U6.36 U1.5
## rowData names(0):
## colnames(6381): AAACCTGAGAAACGCC AAACCTGAGCAATCTC ... TTTGTCATCGCCGTGA
##   TTTGTCATCGTCTGCT
## colData names(0):
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):

While the formatting is different, you can see that the same logic as with the SeuratObject class applies: there are multiple assays and different slots designed to store information related to the single-cell dataset. There is additional information regarding how to operate using this class in the SingleCellExperiment package documentation.

0.5.3 Interoperability between formats

Note that, to allow compatibility between different Bioconductor analysis packages, the Seurat package includes functions to convert to and from SingleCellExperiment objects.

For example, the same SCE object can be generated using the SeuratObject we had already created:

setty_sce <- as.SingleCellExperiment(setty_seurat)

setty_sce

## class: SingleCellExperiment 
## dim: 35102 6381 
## metadata(0):
## assays(2): counts logcounts
## rownames(35102): ENSG00000238009 CICP27 ... ENSG00000278817
##   ENSG00000277196
## rowData names(0):
## colnames(6381): AAACCTGAGAAACGCC AAACCTGAGCAATCTC ... TTTGTCATCGCCGTGA
##   TTTGTCATCGTCTGCT
## colData names(4): orig.ident nCount_RNA nFeature_RNA ident
## reducedDimNames(0):
## mainExpName: RNA
## altExpNames(0):

The function makes a more complete SCE object, including both counts and logcounts assays, as well as the same cell-level metadata that we found in the SeuratObject after creation.

Similarly, a SingleCellExperiment object can be converted to a SeuratObject:

setty_seurat <- as.Seurat(setty_sce)

setty_seurat

## An object of class Seurat 
## 35102 features across 6381 samples within 1 assay 
## Active assay: RNA (35102 features, 0 variable features)

0.5.4 Loading a reduced version of the dataset

In order for this tutorial to run smoothly, we have supplied a reduced version of the Setty et al. dataset featuring a random selection of 3K cells, instead of the ~7K in the original expression matrix. This will prevent memory issues and will allow you to run the analyses in a shorter time span.

The reduced dataset is stored as a h5seurat file, which is a file format designed to store SeuratObjects in an efficient manner. It can be loaded using the LoadH5Seurat() function:

setty_seurat <- LoadH5Seurat("data_preparation/setty_seurat.h5seurat")

## Validating h5Seurat file

## Initializing RNA with data

## Adding counts for RNA

## Adding miscellaneous information for RNA

## Adding command information

## Adding cell-level metadata

## Adding miscellaneous information

## Adding tool-specific results

Note that creating an object with the same name will overwrite the previous SeuratObject. Even though we are doing it on purpose in this case, remember to check for it when you are running your own analyses in R to prevent data loss!

Finally, you may remove unnecessary objects from your R session to free memory before you continue:

rm(setty_sce, setty_rep1_matrix)

0.7.1 Cell quality control

Seurat allows you to easily explore QC metrics and filter cells based on any user-defined criteria. The number of genes and counts per cell, which are standard metrics in the field, are already computed by Seurat and stored in the metadata slot:

head(setty_seurat[[]])

##                  orig.ident nCount_RNA nFeature_RNA
## GACAGAGCAGGAATGC  healthyBM       4160         1855
## GCGGGTTGTCGTTGTA  healthyBM       3662         1991
## AGTGAGGCAGCCTGTG  healthyBM       4418         2179
## AGACGTTGTAGGCATG  healthyBM       4419         2089
## GAGTCCGTCAGTGTTG  healthyBM       6885         2836
## GCAGTTATCCTGTACC  healthyBM       3646         2110

To add to this, we will compute the percentage of reads mapping to mitochondrial genes using the PercentageFeatureSet() function. We will add this information to a new column called percent.mt using the [[ operator, which adds columns to metadata. All mitochondrial gene IDs start with “MT”, which makes it easy to select the genes of interest:

setty_seurat[["percent.mt"]] <- PercentageFeatureSet(setty_seurat, 
                                                     pattern = "^MT-")

We can visualize several QC metrics using the VlnPlot() function:

VlnPlot(setty_seurat, 
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), 
        ncol = 3)

We can see that there are a few outlier cells with high numbers of counts/genes and high mitochondrial gene counts.

To better understand the relationship between the QC metrics, the FeatureScatter() function can be used to visualize feature-feature relationships. We will plot the metrics to find thresholds for data filtering:

FeatureScatter(setty_seurat, feature1 = "nCount_RNA", feature2 = "percent.mt")

FeatureScatter(setty_seurat, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")

Observe the following trends in these plots: - nCount_RNA vs percent.mt: there are a number of cells for which the % of mitochondrial counts is high, whereas the number of total counts is low. These are likely to be dying cells. - nCount_RNA vs nFeature_RNA: some cells show a high number of counts as well as a high number of features (genes). These are likely to be doublets.

Based on the QC plots, we will filter cells with >10% mitochondrial counts and >6000 and <500 genes. This can be achieved using the subset() function in Seurat:

setty_seurat <- subset(setty_seurat, 
                      subset = nFeature_RNA > 500 & 
                      nFeature_RNA < 6000 &
                      percent.mt < 10)

setty_seurat

## An object of class Seurat 
## 35102 features across 2983 samples within 1 assay 
## Active assay: RNA (35102 features, 0 variable features)
##  1 dimensional reduction calculated: pca

As a result, we are left with 2864 high-quality cells for downstream analysis.

0.7.2 Feature quality control

Feature QC is not included in the standard Seurat pipeline -instead, the tool is designed to preserve all genes in the expression matrix, filtering them only upon SeuratObject creation, as was demonstrated before. However, it is often better to to remove genes after cell QC, given that the main filtering criteria is the frequency of non-zero expression across the dataset.

We will calculate the percentage of cells in which is gene is expressed to determine whether the genes should remain in the dataset. To do this, we will need to extract the count matrix from the SeuratObject:

counts <- GetAssayData(setty_seurat,
                       assay = "RNA", slot = "counts")

total_expr <- rowSums(counts > 0)

head(total_expr)

## ENSG00000238009          CICP27 ENSG00000268903 ENSG00000241860 ENSG00000228463 
##               2               3               4               2               4 
## ENSG00000290385 
##               2

Briefly, counts>0 returns a matrix where each entry is TRUE/FALSE if that entry of the counts matrix exceeds 0. Performing rowSums() on that matrix gives the total # of cells in which expression is non-zero for each gene.

We will keep genes that are expressed in at least 3 cells, which constitutes 1% of our current dataset. We will filter the total_expr vector and use the names() function to extract the IDs of our genes of interest:

keep_genes <- names(total_expr[total_expr > 3])

length(keep_genes)

## [1] 25228

head(keep_genes)

## [1] "ENSG00000268903" "ENSG00000228463" "ENSG00000230021" "MTND1P23"       
## [5] "MTND2P28"        "MTCO1P12"

The subset() function we used above can filter genes based on a list of IDs supplied to the features argument:

setty_seurat <- subset(setty_seurat, 
                       features = keep_genes)

setty_seurat

## An object of class Seurat 
## 25228 features across 2983 samples within 1 assay 
## Active assay: RNA (25228 features, 0 variable features)
##  1 dimensional reduction calculated: pca

As a result, 26551 genes remain in the dataset.

0.8.1 Normalization

After removing unwanted cells from the dataset, the next step is to normalize the data. By default, Seurat employs a global-scaling normalization method, “LogNormalize”. This method normalizes the feature expression measurements for each cell by the total expression, multiplies this by a scale factor (10,000 by default), and log-transforms the result. While there are other, more complex normalization strategies, this will suffice to remove library size biases from our data.

Normalized values are stored in setty_seurat[["RNA"]]@data, whereas the raw expression matrix is preserved in the setty_seurat[["RNA"]]@counts slot.

setty_seurat <- NormalizeData(setty_seurat, 
                             normalization.method = "LogNormalize")

setty_seurat[["RNA"]]@data[1:4, 1:4]

## 4 x 4 sparse Matrix of class "dgCMatrix"
##                 GACAGAGCAGGAATGC GCGGGTTGTCGTTGTA AGTGAGGCAGCCTGTG
## ENSG00000268903          .               .                .       
## ENSG00000228463          .               .                .       
## ENSG00000230021          .               .                .       
## MTND1P23                 3.00789         1.866092         2.824224
##                 AGACGTTGTAGGCATG
## ENSG00000268903                .
## ENSG00000228463                .
## ENSG00000230021                .
## MTND1P23                       .

0.8.2 Feature selection

Next, we will calculate a subset of features that exhibit high cell-to-cell variation in the dataset (i.e, they are highly expressed in some cells, and lowly expressed in others). It has been described that focusing on these genes in downstream analysis helps to highlight biological signal in single-cell datasets.

In Seurat, this is done using the FindVariableFeatures() function, which stores results in the setty_seurat[["RNA"]]@var.features slot:

setty_seurat <- FindVariableFeatures(setty_seurat)

# view top 10
head(setty_seurat[["RNA"]]@var.features, n = 10)

##  [1] "HBB"    "CLC"    "LYZ"    "ELANE"  "LTBP1"  "AHSP"   "PRTN3"  "H4C3"  
##  [9] "S100A8" "LTB"

Seurat also includes a function, VariableFeaturePlot() to plot variable features in order to better understand their properties:

VariableFeaturePlot(setty_seurat)

These features will be used in downstream analysis, including PCA. We run a previous PCA with all genes for exploratory purposes, however, it is generally more efficient to select a set of highly variable features, as they capture the majority of data variation. Note that data will need to be re-scaled using the normalized counts instead of the raw counts:

setty_seurat <- ScaleData(setty_seurat, 
                          features = all_genes)

## Centering and scaling data matrix

setty_seurat <- RunPCA(setty_seurat, 
                       features = VariableFeatures(setty_seurat),
                       verbose = FALSE)

DimPlot(setty_seurat)

0.10.1 Determining the dimensionality of the data

To overcome the extensive technical noise in any single feature for scRNA-seq data, Seurat clusters cells based on PCA results. The top principal components represent the majority of the variation in the data, and therefore constitute a robust compression of the dataset. However, how many components should we choose to include?

A common method to select the number of relevant PCs is ranking them based on the percentage of variance explained by each one, which is known as an elbow plot and can be generated by Seurat using the ElbowPlot() function:

ElbowPlot(setty_seurat)

The “elbow” will be the point in the plot where the percentage of explained variance starts to flatten, and in which the majority of the signal in the data will already be captured. In this case, the elbow seems to be situated around 6 PCs, but anything up to 10 PCs might be a good choice, since there is a slight plateau after PC 6. We will select 10 PCs for clustering.

0.10.2 Clustering the cells

To find cell clusters, Seurat first computes the “distances” (i.e. the similarity) among cells in the dataset using the FindNeighbors() function and the number of selected PCs, supplied to the dims argument. Then, groups of cells that lay very close together are detected using the FindClusters() function, which requires users to define a resolution parameter. Higher resolution values will lead to producing a larger number of clusters. For our 3K cell dataset, we will set a resolution value of 0.5:

setty_seurat <- FindNeighbors(setty_seurat, dims = 1:10)

## Computing nearest neighbor graph

## Computing SNN

setty_seurat <- FindClusters(setty_seurat, resolution = 0.5)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2983
## Number of edges: 95351
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8756
## Number of communities: 10
## Elapsed time: 0 seconds

Seurat generated 11 cluters. The clusters can be found using the Idents() function, however, identities are also stored in the seurat_clusters column of the metadata table:

head(Idents(setty_seurat))

## GACAGAGCAGGAATGC GCGGGTTGTCGTTGTA AGTGAGGCAGCCTGTG AGACGTTGTAGGCATG 
##                6                3                0                0 
## GAGTCCGTCAGTGTTG GCAGTTATCCTGTACC 
##                7                2 
## Levels: 0 1 2 3 4 5 6 7 8 9

table(setty_seurat$seurat_clusters)

## 
##   0   1   2   3   4   5   6   7   8   9 
## 592 519 394 344 276 233 220 179 163  63

Let’s see how these clusters are distributed in the PCA plot,

DimPlot(setty_seurat, 
        reduction = "pca",
        group.by = "seurat_clusters",
        label = TRUE)

We observe that we generated several large groups that are clearly different from each other, while some of the smaller clusters are more similar, e.g. clusters 0 and 3. There are also some groups that may not be as clearly defined, such as cluster 8. This type of plot may serve to inform the clustering and re-run it with more appropriate parameters, such as different number of PCs, smaller/larger resolution, etc.

We can also generate a tSNE plot of the clustered cells using the PCs that we used for clustering:

setty_seurat <- RunTSNE(setty_seurat, 
                        dims = 1:10)

DimPlot(setty_seurat, 
        reduction = "tsne",
        label = TRUE)

NOTE: distances between clusters in the tSNE plot do not represent the real similarity between the clusters.

0.10.3 Clustering quality control

Although the groups seem fairly well-defined, it is always interesting to run some quality control analyses to verify that we have correctly clustered the data.

One good way to see whether unwanted factors are influencing the clustering is to check the distribution of QC metrics, such as the number of reads, genes and the proportion of mitochondrial expression. Some of the trends might be biological, but others could be coming from the noise in the data, and it is good to be aware of them in downstream analysis.

VlnPlot(setty_seurat, features = "nCount_RNA")

VlnPlot(setty_seurat, features = "nFeature_RNA")

VlnPlot(setty_seurat, features = "percent.mt")

VlnPlot(setty_seurat, features = "S.Score")

VlnPlot(setty_seurat, features = "G2M.Score")